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This paper discusses stresses and deflections in anisotropic solid struc- 
tures of revolution. It presents two methods based on finite-element tech- 
niques: one, a solid-of-revolution method in which material properties, 
applied forces, and temperatures are independent of angle, and two, a 
long-cylinder method in which these properties are independent of the longi- 
tudinal coordinate. These methods are postulated on uniform stress fields 
within the element, rather than on the usual functional displacement 
description within the element. A Fortran program has been written for 
both these methods, and ample test problems are presented to validate the 
methods. An application is presented for thermal stresses induced during 
the post-growth cooling stage of Czochralski-grown lithium tantalate 
crystals. 

I. INTRODUCTION 

This paper considers two finite-element networks. The first network 
consists of triangular annuli forming a solid of revolution of any 
arbitrary cross section in the radial-longitudinal plane. The network 
has 2t symmetry with regard to material properties, external loadings, 
and temperatures. The second network consists of trapezoidal and 
triangular elements in the plane of the circle forming a right circular 
(actually, a polygon cross section) cylinder long in the longitudinal 
direction. The restriction of 2ir symmetry is lifted for the second 
network and is replaced with the restriction that material properties, 
external loadings, and temperatures are independent of the longitudinal 
direction of the cylinder. 

Both methods were programmed on the IBM 370 computer. For 
the first method, two test plane strain problems are presented: one, a 
hollow cylinder subjected to a negative radial pressure and two, a 
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solid cylinder subjected to a linear radial temperature gradient. Com- 
parisons of stresses and deflections with known plane strain solutions 
are excellent in both cases. In the second method, the solution of a 
long cylinder subjected to a linear radial thermal gradient is presented. 
Comparison with a known theoretical solution again was excellent. 

Results are presented for thermal stresses induced in a lithium tan- 
talate crystal during the post-growth cooling stage. Lithium tantalate 
is in crystal class C 3v . By aligning the trigonal axis of the crystal with 
the longitudinal axis of the cylinder, the problem can be analyzed by 
both methods. Comparison of thermal stresses obtained from the two 
methods was very good. 

II. SOLID-OF-REVOLUTION METHOD 

The basic element for this method is the triangular annulus shown 
in Fig. 1. The element, defined in the radial-longitudinal plane (r — z), 
has 27r symmetry. A network of these elements will comprise any 
desired solid of revolution, whether it is solid or hollow, cylindrical or 
conical, or any combination thereof. Material properties, temperature 
distributions, and external forces must be independent of angle. The 
material properties are then limited to an orthotropic system with 
isotropic properties in the plane of the circle. This limitation is lifted 
in the case of the long cylinder method presented in the next section. 




A = V, L,, x L,, 



V = 27rA(r. + r,+ rJ/3 



|"k» 



e ,i = L ii /L ii 



*.. " e :: - e, b i; = e i; - e h 



Fig. 1— Properties of triangular annulus. 
828 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1974 



Solid triangular elements have been used successfully by other 
authors; 1 - 2 however, a different approach is presented here. Rather 
than to assume a displacement function, 1 - 2 a simpler and more direct 
approach is to postulate uniform stress fields in the annulus. This 
method was applied successfully in the mid-60's by David B. Hall for 
the triangular membrane, but unfortunately his work has not been 
published. Hall's triangular membrane is extended here to a three- 
dimensional model by incorporating a uniform hoop stress in the 
annulus. 

The initial strategy is to relate equilibrium between forces at the 
three points of the triangle i, j, k, and stress fields parallel to the sides 
of the triangle o-.y, *,-*, cr*„ and the hoop stress a e . Before we do this, 
let us first compute a few properties of the triangular annulus. These 
properties include the area of the triangle A, length of each side L, y , 
etc., volume V, areas A„, etc. associated with the stress fields an, 
etc., and direction cosines a tJ , b ij} etc. These properties are presented 
in Fig. 1. 

In matrix notation, equilibrium for annulus n between grid point 
forces and the four stress fields described above is given as follows. 

\Fn\ = [5]{«r }, (1) 

where 



[F n ] = [F ir F it F jr F jz F kr F kt \ y 

[(To) = {(Tij(Tjk(rki<T0}, 



and 



Aijbn 





Akibki 


2ttA/3 


A ij&ij 





Akidki 





Aijbij 


— Ajkbjk 





2tA/3 


Aijdij 


— Ajkdjk 











Ajkbjk 


—Akibki 


2irA/3 





Ajkdjk 


— Akidki 






[*]- 



The r and z subscripts attached to the forces F ir , F it , etc. refer to 
the radial and longitudinal directions, respectively. 

The skew stress field an, a jk , aki is merely an artifice to allow us to 
readily establish the equilibrium relationship. We are actually inter- 
ested in the orthogonal field a z , a r , r Tt , as shown in Fig. 2. The relation- 
ship between the two fields including the hoop stress is given below. 



{*.} - U>1M, 



(2) 



where 



\a t a r T rM ae\ 
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and 



[£] = 



Fig. 2— Stress fields. 



\<$i a% alt 

b% b% bh 

aijbn djkbjk cikibki 

1 



After inverting eq. (2) and substituting it into eq. (1), the following 
relationship is obtained. 

{F n \ = m{v»}, (3) 

where 

[F] = IBTDJ-K 

A conjugate relationship to eq. (3) by an application of the principle 
of virtual work can be stated as 



where the strains are 



{*»} =ym i iU n }, 



{tn} = {e z £ r 7rz«fl}, 



(4) 



the deflections conjugate to {F n } are 

{£/„} = {UirUi.UjrUj.UkrUt.}, 

and the superscript t denotes a transposition. 

The stress-strain relationship, including the thermal strains fadT, 
is given as 

M = [C]({ €n } - {/odrj), (5) 
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where [C] is a symmetric 4X4 stress-strain matrix such that 

Cn = C12, Cu = C24, Cu = C23 = C43 = 0, 
and 

JadTl = ja t dT f a T dT / a e dT ■, a T = at. 

Substituting eq. (5) into eq. (4) obtains the stresses in terms of the 
unknown deflections and known thermal strains as 



I*.} = CC](pOTM^}-|/«drJ) 



(6) 



Substituting eq. (6) into eq. (3) obtains contributions to the network 
stiffness matrix and thermal load vector for annulus n or 

{Fn\ =[Knl\U n } - \E n \, (7) 

where 

is a symmetric 6X6 stiffness matrix and 



{E n } =LFTCl\jadT 



is a 6 X 1 thermal load vector. 

Annuli contributions to the network stiffness matrix and thermal 
load vector are additive, and hence eq. (7) can be written for the net- 
work as 

|F) = ZKHU) - \E), (7a) 

where 

[JC] - L [#»] and \E) = E {E n \, 

n n 

and where the number of equations equals twice the number of points 
of the network (one longitudinal and one radial degree of freedom per 
point). 

The [if] matrix in eq. (7a) is singular, as there is no constraint to 
prevent rigid body motion in the longitudinal direction. Rigid-body 
motion in the radial direction is prevented by the hoop stress field, as 
can be seen in eq. (1). (Note that in Eq. (1) <r e is the only non-self- 
equilibrating stress field.) In addition to at least one reference longi- 
tudinal constraint, radial constraints must be provided for solid 
cylinders at points of zero radius to prevent radial (or hoop) motions 
at these singularities. The degrees of freedom can now be partitioned 
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into an unconstrained set denoted by "a" and a constrained set de- 
noted by "b." This partitioning can be represented schematically as 



[K aa K ab ~\ 
K' ab Kbb . 



Ua 

u b = o 



E b 



(8) 



From the upper equation of (8), the previously unknown deflections 
can now be obtained in terms of known external forces and thermal 

loads, or 

[U a \ = lK aa -]-K\F a \ + [E a \). (9) 

From the lower equation of (8) and with the help of eq. (9), the 
forces of constraint, if desired, can be obtained as 

[F b ] = [K ah yiK a ^{F a \ - \K) - ({E b \ 

- ZK ab y[K aa ^[E a }), (9a) 

where {F* b \ represents external forces applied at the constraint points. 
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a = tan _1 [( r o +r i' sin ~2~^ r 0~ r fi 
e, 3 =[-cos(tt-«)i a , sin ia-H)\] 
e 24 = [-cos(a + H)e a , -sin(a+H)e b ] 
71 =[sin 2 ^(r +r i ) 2 + (r +r l ) 2 ]' / ' 

e r = (cos We a , sinH e b ) 
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A = (r n -r.) (r n +r-| sin^- cosA^- (AREA) 



'0 
2c = (r n -r-) cos 



2 



2 — 2 

f = A/4c 



Fig. 3 — Properties of the trapezoidal element. 
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A0.-= 



A", 



e, 2 = [cos (# -^) e a . sin (fl - i^ ) e b ] 

a 12 = *I2 " 'a b 12 = e 12 e b 

^23 = l _s ' n " *«■ cos " %' 

a 23 = «23". b 23 = i 23 T b 

i 3 ,= [- C os(« + ^-)e a . -sin(fl+M)e b ] 

■ai ""^r^ b 3i = »3r' 5 b 
Fig. 4 — Properties of the triangular element. 

III. LONG-CYLINDER METHOD 

In this section, the 2tt symmetry requirements as to material proper- 
ties, external loadings, and temperature are lifted and a two-dimen- 
sional model is constructed in the plane of the circle. The cylinder is 
considered long in the z- (longitudinal) direction and material proper- 
ties, external loadings, and temperatures are assumed independent of z. 
Shear stresses r yt and t„ are assumed to be zero. The most general 
stress-strain relationship considered is the 4X4 submatrix bounded 
by the dashed lines of eq. (18) (appendix) with the additional proviso 
of no coupling between the stresses <r x , <r v , a x , r xy and the strains y vt} 
y tx . Tentatively, we let e, vanish during the initial phase of the analysis. 
This restriction is subsequently removed in a manner similar to that 
described by Timoshenko 3 for long isotropic cylinders. 

The basic building blocks of the long cylinder is the trapezoidal 
element (Fig. 3) and the isosceles triangular element (Fig. 4, solid 
cylinders only). The trapezoidal element is subjected to three constant 
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stress fields: a radial stress a r , a tangential stress <re, and a shear stress 
r r e- The triangular element is likewise subjected to three uniform 
stress fields, each parallel to a side of the triangle. Equilibrium between 
forces collected at the apices of the elements and three uniform stress 
fields are obtained as shown below (see Figs. 3 and 4 for the properties 
of the elements) . 
Trapezoid : 



(10) 



Triangle : 



[F xX 




fa r 


— cae 


Idu) 


FyX 




fbr 


— cbe 


lb l3 


F x2 




fa r 


cae 


— Ian 


t F+ 




fK 


cbe 


— lb 2i 


}F x3 




-fa r 


ca g 


— lau 


Fy3 




-fbr 


cbe 


-lb l3 


F xi 




— fa r 


— cae 


la-iK 


fy* 




-fbr 


— cbe 


Ian 


F xl ] 




—gau 





gan 


Fyi 




—gbu 





gbn 


F x2 




gai2 


-r B /2 





Fyl 




gbu 


-Ta/2 





F xS 







r a /2 


— ff«31 


_FyZ 







r a /2 


— 0&31 




(11) 



The relationship between the orthogonal stress field (<r r <r«T r <0 and 
the skew stress field (an <r 23 o- 3 i) is given below for the triangle. 




(76 = 




cos 2 A0/2 



cos 2 A0/2 




sin 2 A0/2 1 sin 2 A0/2 

-sin A0/2cos A0/2 sin A0/2 cos A0/2J 




(11a) 



After inverting eq. (11a) and substituting it into eq. (11), equilibrium 
is established between the orthogonal stress field and the forces (de- 
noted as \F n }) at the apices of the element or the points of the network. 
Denote this relationship as: 

[Fn] = LH n l{<r b }. (12) 

Similarly, for the trapezoid, eq. (12) will be the shorthand matrix 
notation for eq. (10). 
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Conjugate to the forces {F n } are deflections { U n ) and therefore the 
conjugate relationship to eq. (12) can be constructed as 

\t h \ = \/AlH n y\U n ), (13) 

Where A is the volume (area times unit thickness) of the trapezoid or 
triangle. 

Upon substituting eq. (13) into eq. (26), the stresses in the element 
in terms of the unknown deflections at the apices and the known thermal 
strains are obtained. 

} - lCl(l/A£H n y[U n ) - lB b -]\jadT\y (14) 

Upon substituting eq. (14) into eq. (12), contributions to the net- 
work stiffness matrix and thermal load vector are obtained for the 
element 

{F n \ =LK n -]\U n } - \E n \, (15) 

where 

[K n -] = l/AlH n XC3lH n y 

is a symmetric 8 X 8 or 6 X 6 stiffness matrix and 

\E n ] = LH^lCTE^ljadT 

is an 8 X 1 or 6 X 1 thermal load vector. 

Contributions of each element to the network stiffness matrix and 
thermal load vector are additive, and hence eq. (7a) can represent the 
force balance at all the points of the long-cylinder network as well as 
those of the solid-of-revolution network. 

Noting the material restrictions outlined in the beginning of this 
section, a network of triangles and trapezoids need only occupy one- 
quarter of a circle with planes of symmetry at 8 = and t/2. Hence, 
tangential deflections must vanish at all points along these planes. The 
previously unknown deflections and constraint forces can now be 
solved in terms of known external forces and thermal loads in the same 
manner as was accomplished in the preceding section [_see eqs. (8) 
(9), and (9a)]. 

We are not quite finished. Recall that we have let e z = throughout 
the cylinder, resulting in an axial stress a z []eq. (23)] applied to the 
ends of the cylinder. If we superimpose a uniform axial stress (a,) 
such that the resulting force on the ends of the cylinder is zero, the 
self-equilibrating distribution remaining on the ends will, by St. 
Venant's principle, 3 give rise only to local effects at the ends. The uni- 
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form axial stress correction to be added to eq. (23) is given as 

[<r s [eq. (23)]dA 



M = - 



I 



(16) 



dA 



This correction results in added radial and tangential strains ob- 
tained from eq. (22). 

MO)) = E u (0)(a,) = [E xz cos 2 + E va sin 2 0)<<r z ) (16a) 

(ee(d)) = Bu(9)W,) = (E Xi sm 2 6 + E y ,cos 2 d)(<x z ). (16b) 



E =10x 10 6 
v =0.25 

a =7.5x 10 _6 /°C 
T = 144R/12 
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-4000 



-8000 



-12.000 



> = CONSTRAINT 

Fig. 5— Thermal stresses plane strain problem. 
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The resulting correction to the radial and tangential deflections are 
next obtained. 

(Ur(r, 0)) = r(e T (6)) = r(E xv cos 2 + E yz sm*d)(a z ) (17) 

(Ue(r, 6)) = [ 6 (r(ee(0)) - (U r (r, 6)))d6 
Jo 

= r J" (E i3 (6) - E 13 (d))(a z )dd 

= r sin5 cos0 (E yz — E xz ). (17a) 

Note that the correction (Ue(r, 6)) vanishes at 6 = and 6 = ir/2, 
agreeing with the stipulated boundary conditions stated previously. 

IV. COMPARISONS WITH THEORETICAL SOLUTIONS 

Two plane strain problems are presented for the solid-of-revolution 
method: one, a linear radial thermal gradient applied to an isotropic 
solid cylinder (Fig. 5) and two, a negative pressure applied to the out- 
side circumference of a hollow cylinder (Fig. 6). Comparisons were 




10PSI 




>= LONGITUDINAL 
CONSTRAINT 



55 85 115 

RADIAL COORDINATE IN INCHES 



AXIS OF CYLINDER 

Fig. 6 — Plane strain solution — thick-walled cylinder. 
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made of the finite element computer results with known theoretical 3 
solutions. The agreement was excellent. 

A thermal problem is presented for the long-cylinder method in 
which an isotropic cylinder is subjected to a constant temperature 
plus a radial thermal gradient (Fig. 7). Comparisons were made of 
the finite-element computer results with a known theoretical solution. 3 
The known solution on page 410 of Ref. 3 was found to be in error. In 
the expression for radial displacement (U), insert (1 — 3^)/(l — v) 
for (1 — 2v) and in the expression for a Z) replace 2r with 2. After the 
above corrections were made to the theoretical solution, excellent com- 
parisons resulted as shown in Fig. 7. 

V. THERMAL STRESSES IN LITHIUM TANTALATE CRYSTALS 

Thermal stresses were computed for a lithium tantalate crystal, 
class C 3 v, during the post-growth cooling stage. The crystal was ana- 
lyzed by the two methods presented in this paper on the IBM 370 
computer. For the solid-of-revolution method (Fig. 8), the model con- 
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Fig. 7 — Thermal stresses — long cylinder. 
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Fig. 8 — Network of triangular annuli — lithium tantalate crystal — solid-of- 
revolution model. 



sists of a right cylinder of radius a = 0.4375 inch and length 6 = 3 
inches plus a cone of length 0.875 inch attached to the far (cold) end. 
Radial constraints are provided at all zero radius points. For the long- 
cylinder method (Fig. 9), the model comprises one-quarter of a circle 
with tangential constraints at all points along the two planes of sym- 
metry (0 = and tt/2). 




0.4375- INCH RADIUS 

10 SPACES 



Fig. 9 — Network of triangles and trapezoids — lithium tantalate crystal — long- 
cylinder model. 
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Elastic properties 4 of lithium tantalate are shown in Fig. 10. Axis 
1 is perpendicular to the vertical mirror plane and axis 3 corresponds 
to the trigonal axis of the crystal. If we strike out the small shear 
coupling terms, Cu, in Fig. 10, then we have met the stress-strain 
relationship [eqs. (5) and (18)] requirements of this paper. The 
coefficients of thermal expansion will be assumed constant in the 1, 2 
directions and piecewise linear along the trigonal axis 3. These co- 
efficients 6 are given below. 



ai = a2 = 21.9 x io-y°c 

a 3 = 8.35 X 10- 6 /°C 
- 3.1 X 10- 8 /°C 


1.16 X 10- 6 /°C 



625°C and higher 

500°C - 625°C 

400°C - 500°C 

room temperature — 400°C. 



The crystal is assumed to be in a stress-free condition at an elevated 
temperature distribution Tl , and elastic material behavior is assumed 
linear from the initial state of strain at temperature distribution Tl 
to the final state of stress and strain at room temperature (26°C). 
This assumption rules out plasticity and stress relaxation. It is felt 
that ignoring these nonlinear effects plus ignoring the initial state of 
stress at Tl does not distort the qualitative picture of the stress 
pattern after the cooling-down process. The elevated temperature 
distribution in degrees centigrade is given below. 

Tl = 1500 - 100 (r/a) 2 - 500 z/b, for solid of revolution 
= 1500 - 100 (r/a) 2 , for the long cylinder. 

The addition of the longitudinal gradient for the solid-of-revolution 
method was found to induce negligible stress in the crystal. 

For the solid-of-revolution method (z-growth), material axes 1, 2, 
and 3 correspond to r, 8, and z, respectively. For the long-cylinder 
method, two cases were investigated: one, z-growth where material 
axes 1, 2, and 3 correspond to x, y, and z, respectively, and two, y- 
growth where material axes 1, 2, and 3 correspond to z, x, and y, 
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Fig 10 — LiTa0 3 stress-strain. 
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Fig. 11 — Thermal stress comparison z-growth LiTaOa crystal. 

respectively. A plot of thermal stresses for the z-growth crystals by 
both methods is compared in Fig. 11. Agreement is very good. A 
comparison of thermal stresses by the long-cylinder method of both 
z-growth and y-growth is shown in Fig. 12 at 6 = 45°. (Typically, 
0-stress variations are less than 5 percent.) The maximum stress 
recorded in absolute value is about 39,000 psi for both z-growth and 
^/-growth. The z-growth (material axis 3 corresponds to x) would yield 
the same results as ^/-growth. This can be inferred from examination of 
the material properties of Fig. 10. 

VI. CONCLUSIONS 

Excellent comparisons were obtained for stresses and deflection in 
isotropic cylinders between the methods described in this paper and 
known theoretical solutions (Figs. 5 to 7). Stress comparisons again 
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Fig. 12 — Thermal stress comparison z-growth and y-growth (at = 45°) LiTaOi 
crystals — long-cylinder method. 

were excellent between the two methods described in this paper for 
z-growth LiTa0 3 crystals (isotropic in the plane of the circle) during 
the post-growth cooling stage (Fig. 11). Comparisons of thermal 
stresses between a y-(or x-) growth LiTa0 3 crystal (both anisotropic in 
the plane of the circle) and a z-growth crystal are presented (Fig. 12). 
Differences in the stress distributions were not great enough to favor 
either growth direction; the more important consideration is to slow 
the cooling process down enough to keep radial thermal gradient to a 
minimum. 
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APPENDIX 

Stress-Strain Relationships — Long Cylinder 

The most general stress-strain relationship considered is the 4X4 
submatrix bounded by the dotted lines of eq. (18) as shown below, 
with no coupling between the stresses (<r x , <t v , <t z , t xv ) and the strains 

(jyz, 7k). 
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(18) 



After inclusion of thermal strains, eq. (18) can be restated as 



(19) 



where 



'(To) = \0 x dy <T Z T x y], 

{to} = {6x ty izyxy), 



UdA = { fa z dT jctydT JadT , 

and [.Co] is the 4X4 submatrix of eq. (18). 

Now consider a rotation about z by the angle 6, and let {<r a \ 
= {<r r <re <Tz T r e] be the radial, tangential, longitudinal, and shear stress, 
respectively. The relationship between this rotated stress field and 
{<T } can be expressed as 



where 



[*] = 



{**} = LB2M, 

cos 2 sin 2 2 sin0 cos0 

sin 2 cos 2 -2 sin0 cos0 

10 

•sin0cos0 sin0 cos0 cos 2 — sin 2 



(20) 



THERMAL STRESS DETERMINATION 843 



Conjugate to eq. (20) is the following relationship between the strains 

{«o} = DB]'[e«}, (20a) 

where 

{e a } = {e r ee e„y r e}. 

Eqs. (19), (20), and (20a) can be combined as 

M = LC a -]{e a } - [£][C ]{adr}, (21) 

where 

After multiplying eq. (21) by [C]" 1 and rearranging, the following 
result is obtained. 

[««} = mi**} + ([5]0- 1 {/^r}, (22) 

where 

Ml = C^a]" 1 = ([^JO-TCo]- 1 ^]- 1 . (22a) 

Tentatively, let t z = 0. This results in a longitudinal stress applied 
at the ends of the cylinder. From the third line of eq. (22), the longi- 
tudinal stress can be obtained as 



<T Z = l/E 33 (E 31 O r + E 32 <T6 + #34 T r9 + j OizdT), 

where, from eq. (22a), 

E n = Ex. cos 2 + E vt sin 2 0, 
E 32 = Ex. sin 2 + E yt cos 2 0, 
E 33 = E S t, 
E iA = 2 sin0 cosfl (E v , - E xt ), 

and E xg , E yt , and E tt are obtained from [Co] -1 - 
From eq. (23) we obtain 



where 
and 



M = U>~]{cx h ) - JOO 1/E„ jaAT 



{<rb\ = {o"r <re T r $} 

1 

rn , ^0 1 o 

L ' 1 -E n /E n -E 32 /E 33 -E 3i /E 33 
1 



(23) 



(24) 
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After substituting eq. (24) into eq. (22) and premultiplying by [D]< 
(remember e, = 0), we obtain 

(■i) = MM + [J^l/adr} , (25) 

where 

[e b ] = {e T €ey r e}, 

UBkl-LDymiDj 

and 

[tf »] = [DJ'P]')" 1 

eos 2 sin 2 —E31/E33 sin0cos0 

sin 2 cos 2 -E32/E33 -sin0cos0 

-2sin0cos0 2sin0eos0 -E 3i /E S3 cos 2 — sin 2 

Let [C] = C^ftU" 1 . After inverting eq. (25) we obtain a desired result. 

{»»} = [C]({e 6 } - [J3 6 ]{ fadT^j (26) 

subject to the restriction e z = 0, which will be removed after an initial 
solution is obtained. 
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